Proyecto Final
Cárdenas Gallardo Paula Daniela y José Manuel Haces López
02-12-2022
Forecasting - Crecimiento de la Economía Mexicana
Proyecto Final | Series de Tiempo
1. Introducción
El objetivo principal del proyecto es realizar forecasts del Producto Interno Bruto de México.
# Librerías necesarias
library(tidyverse)
library(tsibble)
library(feasts)
library(fable)
library(tsibbledata)
library(fpp3)
library(patchwork)
library(timetk)
library(astsa)
library(nlme)
library(plotly)2. Los datos
Se obtuvo el conjunto de datos a partir del archivo
data_gdp.csv, el cual contiene \(114\) registros (filas) y las siguientes
\(9\) variables (columnas):
date\(\rightarrow\) Fecha en cuartos de año, desde 1993 Q1 hasta 2021 Q2gdp\(\rightarrow\) Producto Interno Bruto Real de México, en pesos mexicanos (MXN) (Gross Domestic Product (GDP))consumption\(\rightarrow\) Gasto real de consumo final del sector privado de México, en MXN (Real Private Sector Final Consumption Expenditure)gov_exp\(\rightarrow\) Gasto real en consumo final del gobierno general de México, en MXN (Real General Government Final Consumption Expenditure)imports\(\rightarrow\) Importaciones reales de bienes y servicios (Real Imports of Goods and Services) de México, en MXNexports\(\rightarrow\) Exportaciones reales de bienes y servicios de México, en MXN (Real Exports of Goods and Services)unemployment\(\rightarrow\) Porcentaje de tasa de desempleo de las personas de 15 años y más en México (Unemployment Rate: Aged 15 and Over: All Persons)usd\(\rightarrow\) Promedio de tasas diarias del tipo de cambio de moneda nacional a dólar estadounidense, en MXN (National Currency to US Dollar Exchange Rate: Average of Daily Rates)uncertainty\(\rightarrow\) Índice Mundial de Incertidumbre para México (World Uncertainty Index)
# Importar los datos
data <- read_csv("data_gdp.csv", show_col_type=FALSE) |>
mutate(date=yearquarter(date)) |>
as_tsibble(index=date)
head(data)## # A tsibble: 6 x 9 [1Q]
## date gdp consumption gov_exp imports exports unemploym…¹ usd uncer…²
## <qtr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1993 Q1 2502224. 1523567. 352614. 364603. 355262. 3.5 3.10 0
## 2 1993 Q2 2542759. 1576174. 357229. 393626. 369430. 3.2 3.11 0.0725
## 3 1993 Q3 2516565. 1560920. 344130. 389554 366336. 3.7 3.12 0.260
## 4 1993 Q4 2604024. 1649564. 355999. 429461. 397139. 3.3 3.13 0.291
## 5 1994 Q1 2585847. 1564540. 363414. 432265. 394126. 3.7 3.18 0.0985
## 6 1994 Q2 2693132. 1660768 370829. 461675. 398668. 3.6 3.34 0.263
## # … with abbreviated variable names ¹unemployment, ²uncertainty3. Objetivos
3.1 Investigación
Según el Banco Mundial de México, durante las últimas tres décadas, este país ha tenido un desempeño por debajo de lo esperado en términos de crecimiento, inclusión y reducción de la pobreza en comparación con países similares. La economía tuvo un crecimiento estimado en poco más del \(2.0\)% anual entre 1980 y 2018, lo que limita el progreso en la convergencia en relación con las economías de altos ingresos. La economía mexicana creció \(4.8\)% en 2021, luego de una caída de \(8.1\)% el año anterior debido a la pandemia de COVID-19. En cuanto a al Producto Interno Bruto (PIB), que es la suma del valor (en dinero) de todos los bienes y servicios de uso final que genera un país o entidad federativa durante un período (comúnmente un año o trimestre), se dice que la economía de un país crece cuando su PIB aumenta de un período a otro. Por el contrario, cuando el PIB disminuye se dice que baja la actividad económica. Es sumamente útil para saber si la economía del país está creciendo o no, es decir, si se produjo más o menos que el año anterior. El cambio en el PIB a lo largo del tiempo es uno de los indicadores más importantes del crecimiento económico.
4. Análisis Exploratorio de Datos (EDA)
4.1. Serie de Tiempo por Variable
# Linear plots
# GDP
data |>
ggplot(aes(x=date, y=gdp, color=gdp)) +
geom_line() +
ggtitle("GDP over time") +
ylab('MXN') + xlab('Time') +
theme(
plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic')) 
# consumption
data |>
ggplot(aes(x=date, y=consumption, color=gdp)) +
geom_line() +
ggtitle("Consumption over time") +
ylab('MXN') + xlab('Time') +
theme(
plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic')) 
# gov_exp
data |>
ggplot(aes(x=date, y=gov_exp, color=gdp)) +
geom_line() +
ggtitle("General Government Final Consumption Expenditure over time") +
ylab('MXN') + xlab('Time') +
theme(
plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic')) 
# imports
data |>
ggplot(aes(x=date, y=imports, color=gdp)) +
geom_line() +
ggtitle("Imports over time") +
ylab('MXN') + xlab('Time') +
theme(
plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic')) 
# exports
data |>
ggplot(aes(x=date, y=exports, color=gdp)) +
geom_line() +
ggtitle("Exports over time") +
ylab('MXN') + xlab('Time') +
theme(
plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic')) 
# unemployment
data |>
ggplot(aes(x=date, y=unemployment, color=gdp)) +
geom_line() +
ggtitle("Unemployment over time") +
ylab('%') + xlab('Time') +
theme(
plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic')) 
# usd
data |>
ggplot(aes(x=date, y=usd, color=gdp)) +
geom_line() +
ggtitle("usd over time") +
ylab('MXN') + xlab('Time') +
theme(
plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic')) 
# uncertainty
data |>
ggplot(aes(x=date, y=uncertainty, color=gdp)) +
geom_line() +
ggtitle("Uncertainty over time") +
ylab('Index') + xlab('Time') +
theme(
plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic')) 
4.2. Verificar si hay variables que requieren transformaciones o ajustes:
# consumption
data |>
ggplot(aes(x=date, y=log(consumption), color=gdp)) +
geom_line() +
ggtitle("Log - Consumption over time") +
ylab('MXN') + xlab('Time') +
theme(
plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic')) 
# imports
data |>
ggplot(aes(x=date, y=log(imports), color=gdp)) +
geom_line() +
ggtitle("Log - Imports over time") +
ylab('MXN') + xlab('Time') +
theme(
plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic')) 
# exports
data |>
ggplot(aes(x=date, y=log(exports), color=gdp)) +
geom_line() +
ggtitle("Log - Exports over time") +
ylab('MXN') + xlab('Time') +
theme(
plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic')) 
No hubo mucha diferencia, así que se trabajará con las originales.
# Seasonal plots
# GDP
data |>
gg_season(gdp) |>
ggplotly()# consumption
data |>
gg_season(consumption) |>
ggplotly()# gov_exp
data |>
gg_season(gov_exp) |>
ggplotly()# imports
data |>
gg_season(imports) |>
ggplotly()# exports
data |>
gg_season(exports) |>
ggplotly()# unemployment
data |>
gg_season(unemployment) |>
ggplotly()# usd
data |>
gg_season(usd) |>
ggplotly()# uncertainty
data |>
gg_season(uncertainty) |>
ggplotly()# Boxplot
# GDP
ggplot(data = data, aes(x = quarter(date), y = gdp,
group= quarter(date), color = as.factor(quarter(date)))) +
geom_boxplot() +
ggtitle("GDP per quarter") +
ylab('MXN') + xlab('Quarter') +
theme(plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic'),
legend.title=element_blank()) 
# consumption
ggplot(data = data, aes(x = quarter(date), y = consumption,
group= quarter(date), color = as.factor(quarter(date)))) +
geom_boxplot() +
ggtitle("Consumption per quarter") +
ylab('MXN') + xlab('Quarter') +
theme(plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic'),
legend.title=element_blank()) 
# gov_exp
ggplot(data = data, aes(x = quarter(date), y = gov_exp,
group= quarter(date), color = as.factor(quarter(date)))) +
geom_boxplot() +
ggtitle("Government Final Consumption per quarter") +
ylab('MXN') + xlab('Quarter') +
theme(plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic'),
legend.title=element_blank()) 
# imports
ggplot(data = data, aes(x = quarter(date), y = imports,
group= quarter(date), color = as.factor(quarter(date)))) +
geom_boxplot() +
ggtitle("Imports per quarter") +
ylab('MXN') + xlab('Quarter') +
theme(plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic'),
legend.title=element_blank()) 
# exports
ggplot(data = data, aes(x = quarter(date), y = exports,
group= quarter(date), color = as.factor(quarter(date)))) +
geom_boxplot() +
ggtitle("Exports per quarter") +
ylab('MXN') + xlab('Quarter') +
theme(plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic'),
legend.title=element_blank()) 
# unemployment
ggplot(data = data, aes(x = quarter(date), y = unemployment,
group= quarter(date), color = as.factor(quarter(date)))) +
geom_boxplot() +
ggtitle("Unemployment per quarter") +
ylab('%') + xlab('Quarter') +
theme(plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic'),
legend.title=element_blank()) 
# usd
ggplot(data = data, aes(x = quarter(date), y = usd,
group= quarter(date), color = as.factor(quarter(date)))) +
geom_boxplot() +
ggtitle("MXN change to USD per quarter") +
ylab('MXN') + xlab('Quarter') +
theme(plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic'),
legend.title=element_blank()) 
# uncertainty
ggplot(data = data, aes(x = quarter(date), y = uncertainty,
group= quarter(date), color = as.factor(quarter(date)))) +
geom_boxplot() +
ggtitle("Uncertainty per quarter") +
ylab('Index') + xlab('Quarter') +
theme(plot.title = element_text(color= '#24815C',
size=14,
face='bold'),
axis.title.x = element_text(color='#1B1818',
size=11,
face='bold.italic'),
axis.title.y = element_text(color='#1B1818',
size=11,
face='bold.italic'),
legend.title=element_blank()) 
4.3.Identificar datos atípicos:
# STL decomposition
# GDP
data |>
model(stl = STL(gdp, robust = TRUE)) |>
components() |>
autoplot() |>
ggplotly()# consumption
data |>
model(stl = STL(consumption, robust = TRUE)) |>
components() |>
autoplot() |>
ggplotly()# gov_exp
data |>
model(stl = STL(gov_exp, robust = TRUE)) |>
components() |>
autoplot() |>
ggplotly()# imports
data |>
model(stl = STL(imports, robust = TRUE)) |>
components() |>
autoplot() |>
ggplotly()# exports
data |>
model(stl = STL(exports, robust = TRUE)) |>
components() |>
autoplot() |>
ggplotly()# unemployment
data |>
model(stl = STL(unemployment, robust = TRUE)) |>
components() |>
autoplot() |>
ggplotly()# usd
data |>
model(stl = STL(usd, robust = TRUE)) |>
components() |>
autoplot() |>
ggplotly()# uncertainty
data |>
model(stl = STL(uncertainty, robust = TRUE)) |>
components() |>
autoplot() |>
ggplotly()5. Modelado
- Intentar:
- Transformations/adjustments
- Forecasting using decomposition
- Model combinations
- Handling outliers
- For multivariate models, using dummy variables
- Fourier transformations to account for seasonality
5.1. Train y Test
# Conjunto de Train
train <- data %>%
filter_index(. ~ "2015 Q4")
# Conjunto de Test
test <- data %>%
filter_index("2016 Q1"~.)5.2. Descomposición y ETS
fit1 <- train %>%
model(
Desc_ETS = decomposition_model(
STL(gdp, robust = TRUE),
ETS(season_adjust ~ error("A") + trend("Ad") + season("N"))
)
)
fit1 %>% report()## Series: gdp
## Model: STL decomposition model
## Combination: season_adjust + season_year
##
## ========================================
##
## Series: season_adjust
## Model: ETS(A,Ad,N)
## Smoothing parameters:
## alpha = 0.9998999
## beta = 0.03426769
## phi = 0.9799997
##
## Initial states:
## l[0] b[0]
## 2557866 15668.05
##
## sigma^2: 2768172874
##
## AIC AICc BIC
## 2423.077 2424.065 2438.208
##
## Series: season_year
## Model: SNAIVE
##
## sigma^2: 18394905.6911p <- train %>%
autoplot(gdp) +
geom_line(aes(y = .fitted), data = fit1 %>% augment(), color = "firebrick")
plotly::ggplotly(p)5.3. Mínimos Cuadrados Generalizados
- Este modelo de Mínimos Cuadrados Generalizados con correlaciones Arima se usará para obtener la Significancia Estadística de las variables y descartar las que no lo son.
mod <- gls(gdp ~ ., data=train %>% select(-(date)), correlation=corAR1(form=~1), method="ML")
sum <- summary(mod)
as.data.frame(sum[19])## tTable.Value tTable.Std.Error tTable.t.value tTable.p.value
## (Intercept) 6.605474e+05 1.207815e+05 5.4689472 4.675225e-07
## consumption 5.216109e-01 7.742529e-02 6.7369570 1.977724e-09
## gov_exp 2.400612e+00 3.438593e-01 6.9813802 6.640501e-10
## imports 1.753510e-01 1.147346e-01 1.5283181 1.302358e-01
## exports 5.182668e-01 1.083916e-01 4.7814302 7.447187e-06
## unemployment -1.844479e+04 6.994607e+03 -2.6370014 9.982838e-03
## usd -2.820689e+03 4.724326e+03 -0.5970563 5.520951e-01
## uncertainty 2.979758e+04 2.343101e+04 1.2717155 2.070261e-01
## date 1.694418e+00 1.300246e+01 0.1303152 8.966324e-015.3.1. Modelos Lineales Simples con Errores Arima
data_fit_lm <- train |>
model(
model_consumption = (ARIMA(gdp ~ consumption)),
model_gov_exp = (ARIMA(gdp ~ gov_exp)),
model_exports = (ARIMA(gdp ~ exports)),
model_unemployment = (ARIMA(gdp ~ unemployment))
)
acc <- accuracy(data_fit_lm) |>
select(.model, .type, MAPE, MAE) |>
arrange(MAPE)
acc## # A tibble: 4 × 4
## .model .type MAPE MAE
## <chr> <chr> <dbl> <dbl>
## 1 model_consumption Training 0.674 22511.
## 2 model_gov_exp Training 1.00 33272.
## 3 model_unemployment Training 1.06 34993.
## 4 model_exports Training 1.07 34935.- Con base en la tabla anterior, seleccionaremos los dos modelos con
menor MAPE para predecir con ellos:
- model_consumption
- model_gov_exp
5.4. Regresión Lineal Simple con la variable ‘consumption’
fit2 <- train %>%
model(LM_Consumption = ARIMA(gdp ~ consumption))
p <- train %>%
autoplot(gdp) +
geom_line(aes(y = .fitted), data = fit2 %>% augment(), color = "firebrick")
plotly::ggplotly(p)5.5. Regresión Lineal Simple con la variable ‘gov_exp’
fit3 <- train %>%
model(LM_GovExp = ARIMA(gdp ~ gov_exp))
p <- train %>%
autoplot(gdp) +
geom_line(aes(y = .fitted), data = fit3 %>% augment(), color = "firebrick")
plotly::ggplotly(p)6 Forecasting
Producir predicción un año adelante (hasta 2022 Q2)
6.1. Descomposición y ETS
fc1 <- fit1 %>%
forecast(h = 22)
fc1 %>%
autoplot(test, level = NULL)
6.2. Regresión Lineal con Errores Arima (Consumption)
# Haciendo el test para la regresión
test_lm_consumption <- test %>% select(c(date, consumption))
# HAciendo el Forecast
fc2 <- forecast(fit2, new_data = test_lm_consumption)
# Ploteando la pred
fc2 %>%
autoplot(test, level = NULL)
6.3. Regresión Lineal con Errores Arima (Gov_Exp)
# Haciendo el test para la regresión
test_lm_gov <- test %>% select(c(date, gov_exp))
# HAciendo el Forecast
fc3 <- forecast(fit3, new_data = test_lm_gov)
# Ploteando la pred
fc3 %>%
autoplot(test, level = NULL)
7. Desempeño
Considerar el MAPE - Concatenar todos los resultados en un DF
desemp1 <- fc1 %>%
accuracy(test) %>%
select(c(.model, .type, MAPE, MAE))
desemp2 <- fc2 %>%
accuracy(test) %>%
select(c(.model, .type, MAPE, MAE))
desemp3 <- fc3 %>%
accuracy(test) %>%
select(c(.model, .type, MAPE, MAE))
desemp <- rbind(desemp1, desemp2)
desemp <- rbind(desemp, desemp3)
desemp## # A tibble: 3 × 4
## .model .type MAPE MAE
## <chr> <chr> <dbl> <dbl>
## 1 Desc_ETS Test 3.13 132423.
## 2 LM_Consumption Test 2.39 106839.
## 3 LM_GovExp Test 2.85 119169.8. Referencias Bibliográficas
Banco Mundial. (2022) México: panorama general. Consultado el 24 de noviembre de 2022, de https://www.bancomundial.org/es/country/mexico/overview
INEGI. (2020) Producto Interno Bruto (PIB). Consultado el 24 de noviembre de 2022, de https://cuentame.inegi.org.mx/economia/pib.aspx?tema=e